Trabajo 1 Planificación Territorial

Esta es una guía detallada de como poder realizar el trabajo práctico N°1 del curso Planificación Territorial.

El objetivo es … jardines … ciudad en 15 minutos … ver cuanto se demora a pie.

1 Organización

Importante

Este informe esta pensado en ejecutarse dentro de una carpeta que contenga el trabajo (por ejemplo, trabajo1) y DENTRO DE ESTA una carpeta llamada datos, la cual contendrá todas las capas para el análisis.

A modo de ejemplo, el trabajo debería tener un orden así

Carpeta_del_trabajo/
├── scprit.r
└── datos/
    ├── DEM Los Ríos.tif
    ├── Jardines Junji.shp
    ├── Jardines Integra.shp
    ├── Limite Urbano Censal Valdivia.gpkp
    └── Unidad Vecinal Valdivia.gpkp

En mi caso, la estructura es así

C:/uni/oit/trabajo1/
├── scprit.r
└── datos/
    ├── dem_tif.tif
    ├── layer_jardines_infantiles_junji_2024_20240619120021.shp
    ├── layer_jardines_infantiles_fundacion_integra_20231112011653.shp
    ├── LUC_Valdivia.gpkg
    └── Unidad_vecinal_Valdi.gpkg

Nota

Cabe resaltar que cada usuari@ puede tener el orden que quiera, aquí tendrá este orden para una mejor reproducibilidad.

2 Datos

2.1 Obtención de Datos

Los datos con los que se trabajará serán los siguientes:

• DEM Alos Palsar: El Modelo Digital de Elevación que fue provisto en la clase. Sin embargo, puede descargarse desde IDE Chile.

• Límite Urbano Censal de Valdivia: Provisto en Clase.

• Unidad Vecinal de Valdivia: Provisto en Clase.

• Jardines Infantiles JUNJI: Descargado desde IDE Chile.

• Jardines Infantiles Fundación Integra: Descargado desde IDE Chile.

Advertencia

Aquí debemos evaluar si incluir datos de jardines privados, aunque está abierto a discusión.

2.2 Pre Procesamiento

2.2.1 Librerías

Las librerías serán las mismas que se han visto en las clases. Son las siguientes:

1library(dplyr)
library(ggplot2)
2library(dodgr)
3library(sf)
4library(h3)
5library(terra)
6library(osmdata)
7library(tmap)

1

Paquete para el manejo de datos. Personalmente prefiero cargar tidyverse ya que incluye otras librerías importantes como lubridate (para manejo de fechas) o ggplot2, para visualización.

2

Librería para analizar teoría de grafos. Encontré una guía introductoria (pero en inglés) aquí.

3

Librería que maneja los datos espaciales. Piensen los datos vectoriales como una tabla excel con atributos (como Nombre, Descripción, Region, etc.) y que al final tiene una columna de geometría (usualmente llamada .geo o geometry), la que provee a un dato de una “forma”.

4

Paquete para hacer los hexágonos.

5

Librería para manejar datos ráster.

6

Para obtener datos de elevación y de Open Street Map (OSM).

7

Esta es para hacer cartografía, tal como se hace en QGIS. Pero con esta librería se hacen los mapas interactivos también.

Nota

Se puede hacer click en los números para una “selección” visual.

2.2.2 Situar el Ambiente

Antes de comenzar, dentro del script deberemos setear el ambiente hacia donde están los datos. Revisar Sección 1 para seguir el siguiente ejemplo.

Si el script (scprit.r) se encuentra dentro de la carpeta trabajo1, pero no dentro de la carpeta datos

setwd("datos/")

Si el script no se encuentra dentro de la carpeta trabajo1, se deberá hacer alusión a la ruta completa de los datos. En mi caso:

setwd("C:/uni/oit/t1/datos/")

Si el script se encuentra dentro de la carpeta datos, en teoría no debería hacerse nada.

Advertencia

Windows indica las rutas de las carpetas con el backslash (\) pero R nos los maneja bien. Para indicar la ruta de una carpeta se debe cambiar este símbolo por un slash (/).

2.2.3 Cargar las capas

Una vez seteados en la carpeta correcta, podemos cargar todas las capas y trabajar en base a eso. Para los datos .shp se deberá ocupar st_read() desde el paquete sf y para datos raster la función rast() de terra.

jardinJunji <- st_read("layer_jardines_infantiles_junji_2024_20240619120021.shp")

jardinIntegra <- st_read("layer_jardines_infantiles_fundacion_integra_20231112011653.shp")

dem <- rast("dem_tif.tif")

En el caso de las capas de LUC y Unidad Vecinal de Valdivia, que se encuentran en formato .gpkg, se deberá indicar un argumento extra dentro de st_read(). Esto ya que el formato geopackage es un contenedor que puede almacenar múltiples capas y otros datos relacionados, todo en un único archivo. Es por esto que hay que indicar que capa se quiere acceder. Para ver el listado de capas que tiene un .gpkg, se puede usar la funcion st_layers():

st_layers("LUC_Valdivia.gpkg")

Driver: GPKG 
Available layers:
    layer_name geometry_type features fields crs_name
1 luc_valdivia Multi Polygon        1     12   WGS 84

Podemos ver que el archivo LUC_Valdivia.gpkg tiene solo una capa, llamada luc_valdivia (a través del campo layer_name). En el caso de las unidades vecinales

st_layers("Unidad_vecinal_Valdi.gpkg")

Driver: GPKG 
Available layers:
      layer_name    geometry_type features fields crs_name
1 unidad_vecinal 3D Multi Polygon       57     11   WGS 84

la capa se llama unidad_vecinal. Con esta información podemos cargar correctamente los datos:

lucVald <- st_read("LUC_Valdivia.gpkg",
                   layer = "luc_valdivia")


uvVald <- st_read("Unidad_vecinal_Valdi.gpkg",
                  layer = "unidad_vecinal")

2.2.4 Proyección

Para poder hacer el informe, debemos tener en consideración la proyección de las capas. Para que el flujo funcione, TODAS las capas deben estar en EPSG:4326, el cual corresponder a la proyección mundial WGS 84, el cual se compone de latitud y longitud (en grados).

Datos Vectoriales

Para ver el CRS de un dato vectorial se podría ocupar la función st_crs(), el cual mostrará toda la información de la proyección.

st_crs(uvVald)

Coordinate Reference System:
  User input: WGS 84 
  wkt:
GEOGCRS["WGS 84",
    ENSEMBLE["World Geodetic System 1984 ensemble",
        MEMBER["World Geodetic System 1984 (Transit)"],
        MEMBER["World Geodetic System 1984 (G730)"],
        MEMBER["World Geodetic System 1984 (G873)"],
        MEMBER["World Geodetic System 1984 (G1150)"],
        MEMBER["World Geodetic System 1984 (G1674)"],
        MEMBER["World Geodetic System 1984 (G1762)"],
        MEMBER["World Geodetic System 1984 (G2139)"],
        MEMBER["World Geodetic System 1984 (G2296)"],
        ELLIPSOID["WGS 84",6378137,298.257223563,
            LENGTHUNIT["metre",1]],
        ENSEMBLEACCURACY[2.0]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433]],
    CS[ellipsoidal,2],
        AXIS["geodetic latitude (Lat)",north,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433]],
        AXIS["geodetic longitude (Lon)",east,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433]],
    USAGE[
        SCOPE["Horizontal component of 3D system."],
        AREA["World."],
        BBOX[-90,-180,90,180]],
    ID["EPSG",4326]]

Sin embargo, como esto muestra toda la información, accederemos únicamente al EPSG accediendo con el signo peso ($).

st_crs(uvVald)$epsg

[1] 4326

Esta capa ya se encuentra en el sistema de coordenadas. De querer confirmar esto, se puede comparar con el operador de igualdad (==). Este operador compara si dos valores son iguales. Devuelve TRUE si son iguales y FALSE si son distintos.

st_crs(uvVald)$epsg == 4326

[1] TRUE

Esto confirma que la capa se encuentra en el sistema de proyección correspondiente. Veamos las otras capas

st_crs(lucVald)$epsg == 4326

[1] TRUE

st_crs(jardinJunji)$epsg == 4326

[1] TRUE

st_crs(jardinIntegra)$epsg == 4326

[1] TRUE

Esto indica que todas las capas se encuentran en WGS 84.

Dato Ráster

De la misma manera, terra tiene la función crs() que con el argumento describe = TRUE y accediendo al código (code) con $ podemos comparar de la misma manera

crs(dem, describe = T)$code == 4326

[1] TRUE

Con esto confirmamos que todas las capas están en la misma proyección.

2.2.5 Extensión

El objetivo del trabajo está pensado en las unidades vecinales de Valdivia (uvVald) pero dentro de los límites de la ciudad (lucVald). Visualicemos como se situan todas las capas.

Make this Notebook Trusted to load map: File -> Trust Notebook

A primera vista podemos dos desafíos:

1. La enorme cantidad de establecimientos en jardines infantiles. Esto ocurre ya que ambas campas (Junji e Integra) se encuentran a nivel nacional.

2. También, el LUC nos provee de los límites a trabajar mientras que las unidades vecinales (que son más grandes), proporciona las “particiones” que tendría la ciudad de Valdivia.

Para ambos problemas se pueden ocupar operaciones espaciales, como por ejemplo, el corte o la superposición. El flujo sería:

1. Cortar las Unidades Vecinales con las delimitaciones de Valdivia (dadas por el LUC).

2. A través del LUC de Valdivia, cortar o filtrar ambas capas de jardines infantiles (Junji e Integra).

De este modo, tendríamos todas las capas con las extensiones correspondientes.

2.2.5.1 Recorte

Hay que diferenciar el tipo de operación espacial que se quiere realizar.

Primero que nada, se debe reconocer nuestra capa plantilla (lucVald), la cual delimitará el espacio de interés. Como precaución validaremos que no existan errores topológicos en las capas. Esto puede hacerse con la función st_make_valid() que mantiene polígonos válidos y corrige errores topológicos (como agujeros o superposiciones).

lucVald <- st_make_valid(lucVald)

De igual manera se puede hacer con la otra capa de polígonos

uvVald <- st_make_valid(uvVald)

Ahora, en el caso de los puntos, nos interesa obtener los puntos que se encuentran exactamente dentro del Limite Urbano Censal de Valdivia y a su vez, mantener las propiedades de estos (atributos). En el caso de los puntos esto puede verse así

1junjiVal <- jardinJunji %>%
2  st_filter(lucVald,
3            .predicate = st_covered_by)

1

Crea el objeto junjiVal que es la capa nacional de los jardines JUNJI (jardinJunji). A este último se le realizará una operación (en la próxima línea). Esto puede lograrse con la pipa (%>%).

2

Filtramos (con la función st_filter) la capa original con el LUC (lucVald).

3

El método de filtrado es cubierto por (st_covered_by).

El código anterior se podría traducir a lenguaje humano como:

De la capa de jardines JUNJI de Chile, filtra por los que están cubiertos por el Límite Urbano Censal de Valdivia. El resultado de la operación guardalo en junjiVal.

Aplicamos la misma lógica para la otra capa

integraVal <- jardinIntegra %>%  
  st_filter(lucVald, 
            .predicate = st_covered_by) 

Para manejar los datos de la mejor manera, se unirán ambas capas en un solo objeto, claro que diferenciando entre las fuentes de los datos. Primero veamos los atributos de cada capa

integraVal %>% glimpse()

Rows: 15
Columns: 14
$ FUENTE     <chr> "F. Integra - Agosto 2013", "F. Integra - Agosto 2013", "F.…
$ LAT        <dbl> -39.84270, -39.80373, -39.83497, -39.81447, -39.84400, -39.…
$ LONG       <dbl> -73.21701, -73.20892, -73.25607, -73.22523, -73.22831, -73.…
$ REGION     <chr> "XVI Región de Los Ríos", "XVI Región de Los Ríos", "XVI Re…
$ COMUNA     <chr> "Valdivia", "Valdivia", "Valdivia", "Valdivia", "Valdivia",…
$ NOMBRE     <chr> "160101 El Canelito", "160102 Cuncunita", "160103 Esperanza…
$ MODALIDAD  <chr> "Jardin Infantil", "Jardin Infantil", "Jardin Infantil", "J…
$ DIRECCION  <chr> "Rene Schneider 3210", "Balmaceda 5571", "Calle Arica 2315"…
$ TELEFONO   <dbl> 216531, 224266, 225616, 219868, 232188, 293347, 245301, 273…
$ ZONA       <chr> "Urbana", "Urbana", "Urbana", "Urbana", "Urbana", "Urbana",…
$ MACRONIVEL <chr> "Sala Cuna y Parvulos", "Sala Cuna y Parvulos", "Sala Cuna …
$ CAP_SC     <dbl> 30, 36, 18, 24, 62, 24, 21, 0, 12, 36, 24, 17, 42, 0, 0
$ CAP_PA     <dbl> 174, 71, 96, 96, 118, 15, 64, 62, 29, 32, 32, 0, 32, 0, 0
$ geometry   <POINT [°]> POINT (-73.21701 -39.8427), POINT (-73.20892 -39.80373), PO…

junjiVal %>% glimpse()

Rows: 20
Columns: 14
$ OBJECTID   <dbl> 2731, 2732, 2735, 2736, 2737, 2738, 2740, 2741, 2742, 2743,…
$ N_REGION   <dbl> 14, 14, 14, 14, 14, 14, 14, 14, 14, 14, 14, 14, 14, 14, 14,…
$ REGION     <chr> "LOS RIOS", "LOS RIOS", "LOS RIOS", "LOS RIOS", "LOS RIOS",…
$ COD_COMUN  <dbl> 14101, 14101, 14101, 14101, 14101, 14101, 14101, 14101, 141…
$ COMUNA     <chr> "VALDIVIA", "VALDIVIA", "VALDIVIA", "VALDIVIA", "VALDIVIA",…
$ COD_JUNJI  <dbl> 14101001, 14101002, 14101013, 14101015, 14101016, 14101017,…
$ NOMBRE     <chr> "CAMPANITA", "CARRUSEL", "GOTITAS DE LLUVIA", "FRANCIA", "K…
$ DIRECCION  <chr> "AV. ARGENTINA 2453 POB. VALPARAISO", "AV. DON BOSCO 2701 P…
$ PROGRAMA_E <chr> "JARDÍN INFANTIL CLÁSICO DE ADM. DIRECTA", "JARDÍN INFANTIL…
$ MODALIDAD  <chr> "JARDÍN INFANTIL", "JARDÍN INFANTIL", "JARDÍN INFANTIL", "J…
$ CAPACIDAD_ <dbl> 137, 166, 156, 72, 72, 86, 70, 60, 46, 72, 60, 60, 20, 14, …
$ LONGITUD   <dbl> -73.21223, -73.21348, -73.22359, -73.22923, -73.20537, -73.…
$ LATITUD    <dbl> -39.82491, -39.83058, -39.85030, -39.83890, -39.82980, -39.…
$ geometry   <POINT [°]> POINT (-73.21223 -39.82491), POINT (-73.21348 -39.830…

Ambas capas comparten atributos similares como COMUNA o DIRECCION. Por otro lado, hay atributos que para los objetivos de este trabajo no son necesarios, como latitud, teléfono, etc. Para simplificar esta base de datos compartida nos quedaremos con las siguientes columnas:

• Nombre

• Dirección

• Modalidad

• Fuente: el cual se creará para ver de donde se extrajeron los datos.

Primero se seleccionarán los atributos correspondientes en ambas capas:

1integraValCopia <- integraVal %>%
2  select(nombre = NOMBRE,
3         direccion = DIRECCION,
         modalidad = MODALIDAD) %>%
4  mutate(fuente = "INTEGRA")

junjiValCopia <- junjiVal %>%
  select(nombre = NOMBRE,
         direccion = DIRECCION,
         modalidad = MODALIDAD) %>%
  mutate(fuente = "JUNJI")

1

Aquí, seleccionamos las capa a editar e indicamos que viene una operación con la pipa (%>%).

2

Seleccionamos las columnas de interés con la función (select()). Además, a través de la expresión nombre = NOMBRE se aprovecha de renombrar los nombres originales de las columnas, pasando de NOMBRE a nombre.

3

De la misma forma, se seleccionan las columnas DIRECCION y MODALIDAD además de renombrarlas.

4

Con la función mutate() mutamos la tabla. En este caso, se crea la columna fuente, e indicamos que los valores de esa columna es un texto. Dependiendo de donde sea nuestra fuente, tendrá los valores INTEGRA o JUNJI.

Nota

Es recomendable revisar las capas generadas en cada proceso, ya que este proceso puede ser abstracto y confuso. Puede efectuarse con la función glimpse() o simplemente escribiendo el nombre del objeto. Una forma rápida de hacer esto sería por ejemplo

capa_a_analizar %>% glimpse()

Donde capa_a_analizar es el objeto que se crea, como junjiValCopia o integraValCopia.

Una vez tengamos ambas capas preparadas podemos unir las capas con la función bind_rows(), la cual une tablas por filas. Como requisito, esta función necesita que ambas tablas tengan las mismas columnas (y con los mismos nombres).

jardinesValdivia <- integraValCopia %>% 
  bind_rows(junjiValCopia) 

En este caso, estamos indicando que la capa integraValCopia (que es una tabla), se le una (con bind_rows()) la tabla junjiValCopia. El resultado de esta operación se guarda en jardinesValdivia.

En el caso de las unidades vecinales la opción del filtrado no sirve ya que la capa mantendrá los polígonos que toquen el LUC, pero no los cortará. Esto provocará que la capa de UV sea más grande que el LUC y por ende, no coincidirán espacialmente. Por lo tanto, se debe ocupar la función st_intersection(), la cual corta una capa con otra.

uvValdCortado <- uvVald %>% 
1  st_intersection(lucVald[0]) %>%
2  st_make_valid()

1

Al LUC (lucVald) se le agrega el elemento [0] para deseleccionar todas las columnas (o seleccionar ninguna columna). Así, al revisar uvValdCortado, solo quedarán las columnas del objeto original, es decir uvVald. De no indicar este parámetro, cada polígono de uvValdCortado habría quedado con los atributos de la capa original y de la capa plantilla lucVald.

2

Repetimos la lógica de evitar errores topológicos.

Nota

El nombre de las variables puede ser confuso, como es en el caso de uvValdCortado con uvVald. Estos nombres son arbitrarios y usted puede cambiarlo a su conveniencia. Sin embargo, esta guía seguirá empleando los nombres declarados.

Si desea, puede guardar los objetos creados como un archivo .shp, lo que es recomendable para no tener que repetir el proceso cada vez que se quiera trabajar con estas capas.

st_write(jardinesValdivia, "jardinesValdivia.shp")
st_write(uvValdCortado, "uvValdCortado.shp")
st_write(junjiVal, "junjiVal.shp") # Opcional
st_write(integraVal, "integraVal.shp") # Opcional

Veamos como quedaron las capas con todo el pre proceso.

Make this Notebook Trusted to load map: File -> Trust Notebook